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This is the final report summarizing the work done during the last three years 
under NASA Grant Grant NAG5-8946. Our efforts centered on a systematic development 
of a new generation of three dimensional magneto-hydrodynamic {MHD} numerical 
code, which models the interaction processes of the solar wind or fast flowing 
magnetospheric plasma with “non-magnetic” solar system bodies (e.g. Venus, Mars, 
Europa, Titan). We have also worked on a number of different, more specific and discrete 
studies, as various opportunities arose. In the next few pages we briefly summarize these 
efforts. 


The interaction of the solar wind with “non-magnetic” solar system bodies is very 
different from those that have strong magnetic fields, such as the Earth or Jupiter. The 
obstacle to the supersonic solar wind flow is the ionosphere/atmosphere system at planets 
such as Mars and Venus. A contact discontinuity, called the ionopause, is formed, 
separating the shocked solar wind protons from the planetary thermal ions. Thus the 
nature of the interaction has a strong influence on the structure of the upper ionosphere, 
the formation of the ionopause and magnetic pile-up boundary (also sometimes called 
magnetic barrier). Furthermore the nature, shape and location of the bow shock are 
influenced by the ionosphere, upper atmosphere and exosphere of the planet. In order to 
achieve a thorough understanding/description of the nature and behavior of the plasma 
environment of these solar system bodies, we must be able to elucidate and thoroughly 
understand the physical processes which are responsible for and control these 
interactions. 

The use of gasdynamic and later MHD models to study the interaction of the 
solar wind with the Earth began about three decades ago [e.g. Spreiter et al., 1966]. 
Variations of these models have also been successfully applied to other magnetic planets, 
such as Jupiter [e.g. Ogino et al., 1998; Hansen, 2001]. Both gasdynamic and MHD 
single-fluid models have also been used to study the interaction of the solar wind with the 
“non-magnetic” planets, Venus and Mars [e.g. Spreiter et al., 1970; Murawski and 
Steinolfson, 1996; Bauske et al., 1998]. More recently, a semi-kinetic model has been 
employed to study the solar wind interaction with Venus and Mars [e.g. Brecht, 1997]. 
These semi-kinetic model calculations have been very useful in showing the lack of 
symmetry in the interaction processes. It has been argued that the use of semi-kinetic 
models is especially appropriate for Mars, where the ion gyroradius is of the same order 
as the planetary radius. However, it is important to note that the ideal MHD equations 
have been found to be successful even in such situations. The possible reason for this 
maybe the fact that, as shown by the semi-kinetic code, and observed by the Phobos-2 
wave instrument [Grard et al., 1989], and the MGS magnetometer [e.g. Acuna et al., 
1998; Cloutier et al., 1999], significant wave activity and turbulence are present, leading 
to a wide variety of wave particle interactions, which in turn act as pseudo-collisions. All 
the models mentioned so far have used a solid conducting sphere as the inner boundary in 
their formulation, even though we know that the ionosphere/atmosphere system is the 
actual obstacle. During the last decade multi-species and/or multi-fluid MHD models 
have been introduced to study these interaction processes [Tanaka and Murawski, 1997; 
Liu et al., 1999a; Sauer et al., 1996], providing a significant step forward in these 
calculations. 
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The Aerospace Engineering Department of the University of Michigan has among 
its faculty some of world's best computational fluid dynamicists. A number of years ago 
some of us in the planetary science and space physics community in the Atmospheric, 
Oceanic and Space Science Department recognized the opportunity that the presence of 
these individuals provided and proposed to develop multidimensional MHD codes, with 
an adaptively refined unstructured grid system. Significant funding was received to 
develop such exciting, new numerical codes for space-weather studies. This in turn led us 
to realize that "piggybacking" on these developments we can open the door to a new class 
of MHD models to study Mars, Venus and other non-magnetic solar system bodies at 
significantly reduced additional cost and effort. Thus we began a systematic step by step 
effort, which has lead to the current multi-species, 3D, MHD models, with an ionosphere 
and mass loading processes included. 

During the last grant period we completed the development of the model which 
transformed it from a single to a multi-species formulation. This work was the Ph. D. 
dissertation of one of our graduate students, Ms. Yifan Liu [Liu, 2000], This new 
numerical code is an extensively modified version of a numerical model, named BATS- 
R-US (Block Adaptive-Tree Solar-wind Roe-type Upwind Scheme) [c.f. Powell et al., 
1999], which is the latest version of the single-fluid MHD code developed here at 
Michigan, under separate NASA and NSF sponsorship. The BATS-R-US solution 
method is a highly scalable, massively parallel, block-adaptive mesh refinement (Block- 
AMR) algorithm developed for space physics applications that makes use of recent 
algorithmic advances in high-resolution upwind technology. The basic elements of the 
modified solver are: 1) a cell-centered upwind finite- volume approach that solves the 
hydrodynamic and electromagnetic equations in a tightly coupled manner (rather than in 
separate steps); 2) a flux function based on a new 9-wave approximate Riemann solver 
that accounts for the propagation of MHD waves in a way that is stable, accurate, and 
conservative; 3) limited linear reconstruction that provides second-order accuracy away 
from high-gradient regions and monotonic non-oscillatory solutions throughout the 
computational domain; 4) explicit multi-stage time stepping (source terms are treated 
point-implicitly); 5) block-based solution-adaptive Cartesian grid and data structure; 6) 
physics-based refinement/coarsening; and 7) a parallel implementation that yields 
extremely high computational performance on a variety of massively parallel 
architectures. The resulting parallel adaptive upwind solution algorithm is both robust 
and accurate across a wide range of values of plasma p (the ratio of thermal and magnetic 
pressures). The BATS-R-US numerical scheme has been described in some detail by 
Powell et al., [1999]. Even with the aid of solution-adaptive grids, however, resolving 
three-dimensional, unsteady flows requires massive computing power. This can be 
achieved most efficiently by employing massively parallel computers. The scheme 
incorporated in BATS-R-US is highly scalable, exhibiting excellent parallel performance. 
The code, under the HPCC umbrella, has achieved near perfect scalability and better than 
344 Gflops (this corresponds to better than 20% peak performance) on a 1490-node Cray 
T3E. Good performance has also been obtained on SGI Origin systems and PC Beowulf 
clusters. 
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The new, modified multi-fluid version of the BATS-R-US code solves the MHD 
equations, in a tightly coupled manner. The relevant equations, eigenvectors, propagation 
speeds etc. for the multi-species model have been derived and calculated [Liu., 2000], 
and the model has been coded for multi ion-species implementation. Continuing with our 
philosophy of code development, which is a step-by-step approach, first the simplest two- 
species version was tested and compared to the results from single fluid calculations. This 
multi-species code was then used to investigate the solar wind interaction with Mars [Liu 
et al., 1999a, b, c; Liu, 2000; Liu et al., 2001], Europa [Liu et al., 2000; Liu, 2000] and 
Titan [Liu, 2000; Nagy et al., 2001a], A number of seminars, tutorials and papers were 
also given on the general subject of the interaction of rapidly flowing plasmas with non- 
magnetic solar system objects [e. g. Nagy and Gombosi, 2000; Nagy, 2000a; Nagy and 
Gombosi, 2001; Nagy, 2002] 

In the current 3-species version of the model we solve the following set of MHD 
equations: 
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and where pi, P2 and P3 are the H + , 0 2 + and 0 + mass densities, respectively, S 2 , S 3 and 
L 2 , L3 are the 0 2 + , and 0 + mass source and loss rates, respectively, p is the total thermal 
pressure of the plasma, u is the velocity of the plasma, v is the ion neutral collision 
frequency (taken to be 4xlO‘ 10 { [O] + [C02]} s' 1 , T 0 is the temperature of the newly 
produced ions, 7 is the ratio of specific heats (and taken to be 5/3) and the other symbols 
have their usual definition. 

The Viking 1 and 2 retarding potential analyzer (RPA) measurements [Hanson et 
al., 1977) and subsequent theoretical models [Chen et al., 1978; Fox, 1993] have shown 
that 0 2 + is the major dayside ion below about 300 km in the ionosphere of Mars and 0 + 
becomes important above about 200 km. Thus we selected H + , 0 2 + and 0 + as the three 
ions to use in our three- species model. The chemical reactions that we considered for the 
production and loss of 0 2 + and 0 + ions are: 
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C0 2 + hv — > C0 2 + c 

(7) 

O + hv 0 + + e 

(8) 
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(11) 

0 2 + e — > O + O 

(12) 


We took the photoionization rate of CO 2 and O to be 7.3xl0' 7 and 2.73xl0' 7 sec' 1 , 
respectively [Schunk and Nagy, 2000]; these ionization rates are multiplied by cos(SZA) 
on the dayside, where SZA is the solar zenith angle and no ionization is assumed to be 
present on the nightside. The rates for the reactions shown in equations (9), (10), (1 1) and 
(12) were taken to be 1.64xlO' 10 , 9.6xl0' n , l.lxlO' 9 , and 7.38xl0' 8 cm 3 sec' 1 , 
respectively [Schunk and Nagy, 2000]. We also made the reasonable assumption that 

which simplified the source term. The thermal component of the neutral 

atmosphere was assumed to consist of CO 2 and O; their density distribution was taken to 
be: 

[C0 2 ]= lxl0 10 exp{-(z-140)/ 15.8} and [O] = 3xl0 8 exp{-(z- 140)/ 43.5} 
where the altitudes, z, are in km. 

A computational domain defined by -24R M < x<8R M , - 16R M < y, z<16R M , 
where Rm= 3396 km is the radius of Mars, was used in the calculations and the inner 
boundary was taken to be 140 km above the Martian surface. An adapted grid consisting 
of 17,984 blocks and 1,150,976 computational cells was used in all of the simulations 
with the smallest cells having a dimension of about 53 km near the inner boundary and 
the largest cells, located in the nightside regions far from the planet, having a length of 
about 6792 km. 

We considered three ion species H + , C> 2 + and 0 + , representing the solar wind and 
the two major ionospheric ions, respectively. The inner boundary conditions were 
specified as follows: pi= 0.3p sw , where p sw is the mass density of the undisturbed solar 
wind. This boundary condition is consistent with the expectation of very low H + densities 
in the lower ionosphere, similar to the situation at Venus [e.g. Taylor et al., 1980]. The 
C> 2 + and 0 + densities at the inner boundary were taken to be the photochemical 
equilibrium value. A reflective boundary was used for u; this boundary condition for u, 
results in near zero velocities at the inner boundary and ionospheric velocities of a few 
km s' 1 in the ionosphere, as expected, assuming Venus-like conditions [c.f. Miller and 
Whitten, 1991]. The sum of the electron and ion temperatures at the inner boundary was 
assumed to be 3000 °K and the pressure was set accordingly. The upstream solar wind 
ion and electron temperatures were set to be 5xl0 4 and 3xl0 5 °K respectively. The IMF 
was assumed to be a Parker spiral in the x-y plane with an angle of 56 degrees and a 
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magnitude of 3 nT and the solar wind velocity was selected to be 500 km s' 1 , 
respectively. We picked the solar wind density to be 4 cm' 3 for the nominal case and then 
chose double and half of this value for high and low solar wind pressure simulations, 
respectively. The magnetic field, B, can be set to zero or any arbitrary value at the inner 
boundary; the choice is dictated by the specific problem one wishes to address. Arkani- 
Hamed [2001] developed a harmonic expansion model for the magnetic field of Mars, 
which matches the observations [Acuna et al., 1998]. We were provided with the 
expansion coefficients for a 60-degree expansion [Arkadi-Hamed, private 
communication]. 

The work that we just completed, using the most current version of our model 
addressed three important issues associated with the solar wind interaction with Mars. 

The various results of these studies have been and will be presented at a variety of 
national and international meetings [Ma et al., 2001; Nagy, 2001a; Ma et al., 2002a; Ma 
et al., 2002b; Nagy et al., 2002a] and has been published [Ma et al., 2002c]. 

We first carried out calculations for the nominal solar wind case for which we 
assumed no crustal magnetic field, but included a hot oxygen geocorona. The hot oxygen 
values were taken from the calculations of Kim et al. [1998], which are similar to the 
values obtained by Hodges [2000]. We selected the solar maximum values, in order to 
maximize the possible impact of such a hot atom corona. We assumed that the hot 
oxygen densities are altitude dependent, but do not change with solar zenith angle on the 
dayside and that they are zero on the nightside. Given that we were specifically interested 
in maximizing the dayside results, this assumption is reasonable. We found, within the 
spatial resolution that our current model allows, no meaningful changes in the dayside 
bowshock or ionopause locations, compared to the case with no hot oxygen present in the 
model. We did not include thermal or hot hydrogen in our calculations, but given the 
mass ratios and the relative densities [e.g. Nagy and Cravens, 1988], they are not 
believed to play an important role. In the model we consider photoionization only. It was 
shown by Zhang et al. [1993] that electron impact ionization may at times exceed 
photoionization by up to a factor of four, although Krimsky and Breus [1996] have 
criticized this conclusion. As we were trying to test for the “maximum” effect, in a 
further run we multiplied the subsolar solar cycle maximum photoionization rate of the 
hot oxygen by a factor of four and kept it constant (solar zenith angle independent) over 
the dayside. However, we still could not see any significant changes in the dayside 
bowshock or ionopause locations. The calculations indicated that the presence of a hot 
oxygen corona does not, within the resolution and accuracy of the model, lead to any 
significant effect on the dayside bowshock and ionopause positions, but does result in a 
30% and 15% increase in the trans-terminator and escape fluxes, respectively. 

Next the trans-terminator fluxes and escape fluxes down the tail were calculated. 
First we made these calculations neglecting the effects of the crustal magnetic field, 
considering nominal, high and low solar wind pressure cases. We also calculated these 
fluxes for the nominal solar wind case and included the effect of crustal fields. The 
calculated fluxes for this latter case are given in Table 1. 
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Terminator 
Flux (s' 1 ) 

Escape 
Flux (s ') 

0 2 + 

2.72xl0 25 

1.88xl0 25 

0 + 

0.47xl0 25 

0.36xl0 25 

Total 

3.19xl0 25 

2.24xl0 25 


Table 1. Calculated trans-terminator and escape fluxes for nominal solar wind 
parameters, and the remnant crustal fields considered. 


The calculated flux values were found to be consistent with the escape fluxes of a few 

25 -1 

times 10 sec , estimated from the Phobos-2 measurements [Lundin et al., 1989; 
Rosenbauer et al., 1989] and within the calculated limiting fluxes from the dayside 
ionosphere of (4.4-14. 2)xl0 25 sec' 1 [Fox, 1997]. 

Finally, we ran our numerical model to assess the effects of the crustal magnetic 
fields, using the nominal solar wind and atmospheric parameters, with the 60-order 
harmonic field model [Arkani-Hamed, 2001]. The calculations were carried out for the 
case where the location of the largest crustal field (180 longitude) was pointing sunward, 
and the solar wind flow direction was parallel with the equator of the planet. Figure 1 
shows the calculated magnetic field and velocity components in the equatorial plane, 
while Figure 2 shows the same parameters in the noon-midnight meridianal plane. The 
color plots show the magnitude of the magnetic field and velocity components; the white 
lines show the field line directions and the arrows indicate the direction (not the 
magnitude) of the velocities. The draping of the IMF is clearly visible in the equatorial 
plane as indicated in Figure 1 ; the small “bends and kinks” are the result of the crustal 
fields. There also appears to be a flux rope like feature; flux ropes have been observed by 
the MGS magnetometer [Vignes et al., 2002]. The general trend of the plasma flow is 
upward and toward the terminator on the dayside and either downward into the 
ionosphere (helping to maintain the nightside ionosphere) or down the tail (and escaping 
from the planet) on the nightside. Figure 2 shows the presence of closed field lines (mini 
“magnetocylinders”), the result of the “merging” of the crustal and IMF fields. The 
presence of such mini-magnetospheres have been inferred by the electron reflectometer 
carried aboard the Mars Global Surveyor (MGS) spacecraft [Mitchell et al., 2001]. 
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Figure 1. The component of the calculated magnetic field and velocity in the x-y plane. 
The color plots show the magnitudes; the white lines marked with arrows 
indicate the vector directions of the magnetic field and the arrows the direction 
(not the magnitude) of the velocity. 
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Figure 2. The component of the calculated magnetic field and velocity in the x-z plane. 
The color plots show the magnitudes; the white lines marked with arrows 
indicate the vector directions of the magnetic field and the arrows the direction 
(not the magnitude) of the velocity. 




10 


During the grant period a number of "studies of opportunities" were also carried 
out with the support of the current grant. We have worked on a variety of calculations 
concerning hot atom populations at Mars and Venus. The most recent effort, carried out 
during the current grant period, was aimed at calculating the hot carbon densities at Mars. 
The calculation used a two-step process: first a two-stream transport code was used to 
solve for the distribution function at the exobase and then these results were used in a 
Liouville equation solution above the exobase. It was found that photodissociation of 
carbon monoxide is the largest source of hot carbon atoms in the upper atmosphere of 
Mars, larger than dissociative recombination of CO + and much larger than the creation of 
hot carbon through collisions with hot oxygen atoms. These conclusions indicate that the 
situation at Mars is different from that at Venus, where the latter two processes are 
believed to be the dominant ones. It was also found that the high solar activity carbon 

densities are about an order of magnitude larger than those for the low solar activity case. 

3 23 

The calculated values of 4.7x10 and 3.6x10 cm' near the exobase, for the high and 
low solar activity conditions, respectively, indicate that atomic carbon remains a minor 
constituent throughout the exosphere, compared to hydrogen and oxygen. These results 
were presented at national [Nagy et al., 2000] and international meetings [Nagy, 2001b] 
and published in the Journal of Geophysical Research [Nagy et al., 2001b]. 

We were also involved in the interpretation of the Galileo radio occultation results 
at Callisto. Observations of a significant ionosphere {- 10 4 cm' 3 } at Callisto, under 
certain circumstances, prompted us to evaluate its implication in terms of the presence of 
an atmosphere/exosphere. It was found that a detectable ionosphere was only present at 
the observed location when the trailing hemisphere of Callisto, which is the one that is 
impacted by the corotating plasma of Jupiter’s magnetosphere, was illuminated by the 
Sun. We considered both photo and electron impact ionization, and used a couple of 
different approaches to estimate the surface density. It was rather surprising and 
reassuring to find that all of the methods used to estimate the surface neutral density gave 
very similar results; the estimated values fall between 1 and 3 x 10 10 cm' 3 . These and 
related radio occultation results were presented at two meetings [Kliore et al., 2001a; 
Kliore et al., 2001b] and a paper is in press in the Journal of Geophysical Research 
[Kliore et al., 2002 ]. 

We were one of the three organizers of the Yosemite Meeting on “Comparative 
Solar System Aeronomy”. The purpose of this meeting was to bring together the various 
communities involved in similar studies of the upper atmospheres and ionospheres of 
different solar system bodies. This attempt at cross fertilization of scientists, often 
working on similar problems and yet not aware of each other’s work was highly 
successful. It brought together about 70 scientists and a number of graduate students with 
the financial assistance of both NASA and NSF. The meeting lasted 5 days and allowed a 
wide range of topics to be covered, including a presentation on solar system ionospheres 
by the P.I. of this proposal [Nagy, 2000b]. This meeting provided plenty of opportunities 
for timely interchanges among the attendants. A book, containing most of the material 
covered at the meeting, as well as additional related topics, has now been published 
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[Mendillo, Nagy and Waite, 2002]; this includes a chapter on solar system ionospheres 
co-authored by the P.I. of this proposal [Nagy and Cravens, 2002]. 

The P.I. of this proposal also participated in the international meeting “Interaction 
of the Solar Wind with Mars; Comparison between MGS and Phobos Results”, organized 
and held at ISSI in Bern, Switzerland, October 22-26, 2001. The purpose of this meeting 
was to bring together an international community of scientists who have worked on the 
Phobos and MGS results, along with modelers. This was a 5-day meeting providing a 
platform for formal presentations [e.g. Nagy, 2001a] and wide ranging discussions. The 
meeting is being followed by the preparation of three extensive review papers and a 
number of contributed ones, to be published in Space Science Reviews and a special ISSI 
book. One of the review papers is “assembled” by the P.I. of this proposal [Nagy et al., 
2003]. 


In summary we supported 2 graduate students and will have presented 
approximately 15 papers at national and international meetings and published 8 papers in 
refereed journals with the full or partial support of the current grant, during its three year 
period. 
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